Elranatamab & Teclistamab IPTW Analysis

1 SETUP

1.1 Open libraries

library(gtsummary)
library(dplyr)
library(tidyverse)
library(survival)
library(survminer)
library(lubridate)
library(readxl)
library(janitor)
library(ggplot2)
library(tidycmprsk)
library(ggsurvfit)
library(tidymodels) 
library(sjPlot)
library(patchwork)
library(bstfun)
library(ComplexUpset)

library(writexl)
library(rms)

library(timeROC)
library(rms)
library(pec)
library(prodlim)
library(riskRegression)

library(plotly)
library(ggtext)

library(WeightIt)
library(cobalt)
library(survey)

2 DATA WRANGLING

2.1 Read in MM datasets

db_Kansas = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 2.11.2025 - due 5-1-25 (de-identified).xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  filter(!is.na(age_at_bs_ab_initiation)) %>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_Kansas, "Kansas_elra_data.xlsx")

db_CC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/CC bsAb Consortium Submission CCF 6-11-25.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_CC, "CC_elra_data.xlsx")

db_HUMC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/BsAb Consortium Data HUMC.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_HUMC, "HUMC_elra_data.xlsx")

db_MCC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data MCC-deidentified.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )


db_Huntsman = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/Utah_Huntsman_bsAb Consortium Data FINAL_5.4.25.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc)/1000,
    baseline_alc = as.double(baseline_alc)/1000,
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_Huntsman, "Huntsman_elra_data.xlsx")


db_FHCC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/FHCC Bispecifics (5.5.2025).xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    date_of_infection_80 = ymd(date_of_infection_80),
    baseline_anc = as.double(baseline_anc),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )


db_UTSW = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 2.11.2025 1.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3")) %>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles)
  )

# Write to Excel
#write_xlsx(db_UTSW, "UTSW_elra_data.xlsx")


db_MDACC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/MDACC BITE ASH 2025.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1),
    total_number_paraskelatal_plasmacytomas = as.double(total_number_paraskelatal_plasmacytomas)
  )

# Write to Excel
#write_xlsx(db_MDACC, "MDACC_elra_data.xlsx")


db_Roswell = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 5.8.2025_NN.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1),
    total_number_paraskelatal_plasmacytomas = as.double(total_number_paraskelatal_plasmacytomas)
  )

# Write to Excel
#write_xlsx(db_Roswell, "Roswell_elra_data.xlsx")

# Combine all into a single dataframe
combined_db <- bind_rows(
  list(
    HUMC = db_HUMC,
    MCC = db_MCC,
    Huntsman = db_Huntsman,
    FHCC = db_FHCC,
    UTSW = db_UTSW,
    MDACC = db_MDACC,
    Roswell = db_Roswell,
    CC = db_CC,
    Kansas = db_Kansas
  ),
  .id = "site"
) %>%
  mutate(
    prior_bcma_directed_therapy_yes_no = ifelse(!is.na(most_recent_bcma_agent_name),1,0),
    bs_ab_name = case_when(
      bs_ab_name == 1 ~ "tec",
      bs_ab_name == 3 ~ "elra",
      .default = NULL
    ),
    max_crs_grade_0_5 = ifelse(is.na(max_crs_grade_0_5), 0,max_crs_grade_0_5), ## Set NAs to 0
    max_icans_grade = ifelse(is.na(max_icans_grade), 0, max_icans_grade ) ## Set NAs to 0
  )

# Convert all datetime columns to Date
combined_db_clean <- combined_db %>%
  mutate(across(where(lubridate::is.POSIXt), as.Date))

# Write to Excel
write_xlsx(combined_db_clean, "elra_tec_data.xlsx")


# Write to Excel
# write_xlsx(combined_db_clean %>% filter(prior_bcma_directed_therapy_yes_no != prior_bcma_directed_therapy), "elra_data_discordant.xlsx")

## Several institutions do not have any cases of elranatamab

# db_UABMC = read_excel(path = "Data/bsAb Consortium Data 5.8.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_Stanford = read_excel(path = "Data/De-identified Stanford bsAb Consortium Data Template 4.30.2025.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_LCI = read_excel(path = "Data/bsAb Consortium Data LCI.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_COH = read_excel(path = "Data/bsAb Consortium Data Template 2.11.2025_COHv2.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_MUSC = read_excel(path = "Data/USMM bsAb Consortium Data Template 2.11.2025.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_Mayo = read_excel(path = "Data/Bispecific_consortium_20250521 NoPHI.xlsx", sheet = "bsAb_DATA") %>%
#   filter(`bsAb name` == 3)%>%
#   clean_names()

2.2 Calculate CAR-HEMATOTOX score and risk level

combined_db <- combined_db %>% 
  mutate(
    # baseline_anc = as.numeric(baseline_anc),
    # baseline_hgb,
    # baseline_pl_ts,
    # baseline_crp_prior_to_bs_ab_initiation,
    # baseline_ferritin_prior_to_bs_ab_initiation,
    HT_Score = 
      (baseline_anc <= 1.2) +
      (baseline_hgb <= 9.0) +
      (baseline_pl_ts >= 76 & baseline_pl_ts <= 175) +
      (baseline_crp_prior_to_bs_ab_initiation >= 3.0) +
      (baseline_ferritin_prior_to_bs_ab_initiation >= 650 & baseline_ferritin_prior_to_bs_ab_initiation <= 2000) +
      (baseline_pl_ts <= 75) + ## Add an extra point for Plt <=75
      (baseline_ferritin_prior_to_bs_ab_initiation > 2000), ## Add an extra point for ferritin >2000
    
    # Classify the risk based on the score
    HT_Risk = factor(
      case_when(
        HT_Score >= 2 ~ "HThigh",
        baseline_pl_ts <=75 ~ "HThigh", ## Equivalent to at least 2 points
        baseline_ferritin_prior_to_bs_ab_initiation > 2000 ~ "HThigh", ## Equivalent to at least 2 points
        (baseline_anc <= 1.2) + (baseline_hgb <= 9.0) + (baseline_pl_ts >= 76 & baseline_pl_ts <= 175) >=2 ~ "HThigh",
        HT_Score <2 ~ "HTlow",
        .default = "NA"
      ),
      levels = c("HTlow", "HThigh")
    )
)
write_xlsx(combined_db %>% select(site, patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate, baseline_anc, baseline_hgb, baseline_pl_ts, baseline_crp_prior_to_bs_ab_initiation, baseline_ferritin_prior_to_bs_ab_initiation, HT_Score, HT_Risk), "HT_score.xlsx")

2.3 Identify severe infections

combined_db <- combined_db %>% 
  mutate(
    severe_infection = ifelse(infection_severity_73 %in% 1 | infection_severity_76 %in% 1 | infection_severity_79 %in% 1| infection_severity_82 %in% 1,1,0),
    first_severe_infection_date = case_when(
      infection_severity_73 %in% 1 ~ date_of_infection_74,
      infection_severity_76 %in% 1 ~ date_of_infection_77,
      infection_severity_79 %in% 1 ~ date_of_infection_80,
      infection_severity_82 %in% 1 ~ date_of_infection_83,
      .default = NULL
    ),
    first_severe_infection_type = case_when(
      infection_severity_73 %in% 1 ~ number_infection_1_bacterial_viral_fungal,
      infection_severity_76 %in% 1 ~ number_infection_2_bacterial_viral_fungal,
      infection_severity_79 %in% 1 ~ number_infection_3_bacterial_viral_fungal,
      infection_severity_82 %in% 1 ~ number_infection_4_bacterial_viral_fungal,
      .default = NULL
    )
    
    
  )

3 PROPENSITY SCORES

3.1 Regression analysis for treatment allocation

theme_gtsummary_compact()

## Setting theme "Compact"

variables <- c("elra", "Age", "Female sex", "ECOG 2+", "EMD", "Hgb", "LDH", "Prior BCMA")

data1 <- combined_db %>% 
  drop_na(age_at_bs_ab_initiation, ecog_ps_at_bs_ab_initiation_0_5, just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk, baseline_hgb, baseline_ldh_prior_to_bs_ab_initiation) %>%
  transmute(
    Age = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female,
    elra = ifelse(bs_ab_name == "elra", 1, 0),
    "ECOG 2+" = ifelse(ecog_ps_at_bs_ab_initiation_0_5 >=2,1,0),
    EMD = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,1,0),
    Hgb = as.double(baseline_hgb),
    LDH = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no
  )

uv_tab <- tbl_uvregression(
    data1[c(variables)],
    method = glm,
    y = "elra",
    method.args = list(family = binomial),
    exponentiate = TRUE
  ) %>%
  sort_p()


mv_tab<-glm(elra ~ Age + `Female sex` + `ECOG 2+` + EMD + Hgb + LDH + `Prior BCMA`, data=data1, family=binomial) %>%
  tbl_regression(exponentiate=TRUE)

tbl_merge(list(uv_tab, mv_tab), tab_spanner = c("**Univariate**", "**Multivariable**"))

Characteristic	Univariate				Multivariable
	N	OR1	95% CI1	p-value	OR1	95% CI1	p-value
Hgb	371	1.12	1.01, 1.25	0.029	1.14	1.02, 1.28	0.026
EMD	371	1.69	0.95, 2.96	0.069	1.71	0.96, 3.04	0.068
Age	371	1.02	1.00, 1.04	0.087	1.01	0.99, 1.04	0.3
Female sex	371	1.43	0.93, 2.21	0.11	1.43	0.92, 2.24	0.11
Prior BCMA	371	0.76	0.49, 1.17	0.2	0.78	0.50, 1.23	0.3
ECOG 2+	371	1.19	0.75, 1.88	0.4	1.27	0.78, 2.07	0.3
LDH	371	1.00	1.00, 1.00	0.8	1.00	1.00, 1.00	>0.9
1 OR = Odds Ratio, CI = Confidence Interval

3.2 IPTW

IPTW.data <- combined_db %>% 
  drop_na(age_at_bs_ab_initiation, ecog_ps_at_bs_ab_initiation_0_5, just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk, baseline_hgb, baseline_ldh_prior_to_bs_ab_initiation) %>%
  mutate(
    Age = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female,
    bsAb = bs_ab_name,
    "Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
    "ECOG" = case_when(
      ecog_ps_at_bs_ab_initiation_0_5 %in% c("0","1") ~"0-1",
      ecog_ps_at_bs_ab_initiation_0_5 ==2 ~"2",
      ecog_ps_at_bs_ab_initiation_0_5 >=3 ~"3+",
      .default = NULL),
    EMD = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,1,0),
    Hgb = as.double(baseline_hgb),
    LDH = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no
  )

## Use the WeightIt function to generate propensity scores and weights
W.out <- weightit(bsAb ~ Age + `ECOG` + EMD + Hgb + LDH + `Prior BCMA`,
                  data=IPTW.data, 
                  method = "glm",    ## Methods: glm, gbm, cbps, npcbps, ebal, optweight, super, bart, energy
                  estimand = "ATE",
                  stabilize = TRUE)

bal.tab(W.out, stats = c("m", "v"), thresholds = c(m=0.25), binary = "std")

## Balance Measures
##                Type Diff.Adj     M.Threshold V.Ratio.Adj
## prop.score Distance   0.0050 Balanced, <0.25      0.9733
## Age         Contin.   0.0090 Balanced, <0.25      0.9595
## ECOG_0-1     Binary  -0.0121 Balanced, <0.25           .
## ECOG_2       Binary   0.0152 Balanced, <0.25           .
## ECOG_3+      Binary  -0.0027 Balanced, <0.25           .
## EMD          Binary  -0.0034 Balanced, <0.25           .
## Hgb         Contin.  -0.0186 Balanced, <0.25      1.0974
## LDH         Contin.  -0.0276 Balanced, <0.25      1.0961
## Prior BCMA   Binary  -0.0119 Balanced, <0.25           .
## 
## Balance tally for mean differences
##                     count
## Balanced, <0.25         9
## Not Balanced, >0.25     0
## 
## Variable with the greatest mean difference
##  Variable Diff.Adj     M.Threshold
##       LDH  -0.0276 Balanced, <0.25
## 
## Effective sample sizes
##              elra    tec
## Unadjusted 123.   248.  
## Adjusted   113.38 242.39

summary(W.out)

##                   Summary of weights
## 
## - Weight ranges:
## 
##         Min                                  Max
## elra 0.4745 |---------------------------| 2.3134
## tec  0.7811      |-------------|          1.7418
## 
## - Units with the 5 most extreme weights by group:
##                                         
##          123    102     81     91      9
##  elra 1.5458 1.5896  1.653 1.7607 2.3134
##          169     71    274    258    271
##   tec 1.4257 1.4453 1.5215 1.6783 1.7418
## 
## - Weight statistics:
## 
##      Coef of Var   MAD Entropy # Zeros
## elra       0.292 0.222   0.040       0
## tec        0.152 0.115   0.011       0
## 
## - Effective Sample Sizes:
## 
##              elra    tec
## Unweighted 123.   248.  
## Weighted   113.38 242.39

love_plot <- love.plot(W.out,
          binary = "std",
          line = TRUE, 
          abs = TRUE,
          thresholds = c(m=0.2),
          sample.names = c("Unweighted", "IPTW"),
          title = "A)",
          var.order = "unadjusted",
            var.names = c(
              prop.score = "Distance"
          ))+ theme(plot.title = element_text(hjust = 0))
love_plot

weight_df = data.frame(W.out$weights, W.out$treat)

IPTW.data$weights <- W.out$weights
  
IPTW.data$ps <- W.out$ps


write.csv(data1, "IPTW.dataset.csv")

3.3 Scatterplot of weights

Weights_plot <- ggplot(
  IPTW.data, aes(x=ps, y = weights, color = bsAb)) + 
  geom_point(shape = 16) +
  theme_classic() + 
  scale_y_continuous(breaks = seq(0,4 , by = 0.5)) +
  labs(x = "Propensity score", y = "Weights") + 
  ggtitle("B)")
Weights_plot

4 TABLES

4.1 Site contributions

theme_gtsummary_compact()

combined_db %>%
  transmute(
    "Site" = site,
    bs_ab_name,
  ) %>%
  tbl_summary(
    by = bs_ab_name,
    missing = "no",
    sort = list(all_categorical() ~ "frequency"),
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)"
    )
  ) %>%
  bold_labels()

Characteristic	elra N = 1301	tec N = 2771
Site
MCC	53 (41%)	72 (26%)
MDACC	12 (9.2%)	57 (21%)
CC	20 (15%)	47 (17%)
UTSW	10 (7.7%)	40 (14%)
FHCC	7 (5.4%)	24 (8.7%)
Huntsman	9 (6.9%)	22 (7.9%)
Kansas	4 (3.1%)	13 (4.7%)
HUMC	10 (7.7%)	1 (0.4%)
Roswell	5 (3.8%)	1 (0.4%)
1 n (%)

4.2 Baseline characteristics

theme_gtsummary_compact()

combined_db %>%
  transmute(
    #"Site" = site,
    "Age (years)" = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female == 1,
    "Race" = factor(case_when(
      is.na(race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5) ~ NA,
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 0 ~ "White",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 1 ~ "Black",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 2 ~ "AAPI",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 3 ~ "American Indian",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 4 ~ "Other",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 5 ~ NA
    ), levels = c("White","Black","AAPI","American Indian","Other",NA)),
    #ecog_ps_at_bs_ab_initiation_0_5,
    "ECOG" = case_when(
      ecog_ps_at_bs_ab_initiation_0_5 == 0 ~ "0",
      ecog_ps_at_bs_ab_initiation_0_5 == 1 ~ "1",
      ecog_ps_at_bs_ab_initiation_0_5 == 2 ~ "2",
      ecog_ps_at_bs_ab_initiation_0_5 >= 3 ~ "3+"
    ),
    "LDH (U/L)" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Hgb (g/dL)" = baseline_hgb,
    "Heavy chain" = factor(case_when(
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgG" ~ "IgG",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgM" ~ "IgM",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgA" ~ "IgA",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgD" ~ "IgD",
      .default = "None"
      ), levels = c("IgG", "IgA","IgM","IgD","None")),
    "True EMD" = just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,
    "Oligo or non-secretory" = oligo_or_non_secretory_at_bs_ab_initiation_0_no_1_yes,
    "≥1 HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes==1 | t_4_14_prior_to_bs_ab_0_no_1_yes == 1 | t_14_16_prior_to_bs_ab_0_no_1_yes == 1 | amp_1_q_only_0_no_1_yes == 1 | gain_amp1q21_prior_to_bs_ab_0_no_1_yes == 1,1,0),
    "del(17p)" = deletion_17p_prior_to_bs_ab_0_no_1_yes,
    "t(4;14)" = t_4_14_prior_to_bs_ab_0_no_1_yes,
    "t(14;16)" = t_14_16_prior_to_bs_ab_0_no_1_yes,
    "gain/amp(1q)" = gain_amp1q21_prior_to_bs_ab_0_no_1_yes,
    "amp(1q)" = amp_1_q_only_0_no_1_yes,
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "Prior ASCT" = prior_auto_sct_no_0_yes_1,
    "Prior CAR-T" = prior_cart_0_no_1_yes,
    "Bortezomib refractory" = bortezomib_status == 1,
    "Lenalidomide refractory" = lenalidomide_status == 1,
    "Pomalidomide refractory" = pomalidomide_status == 1,
    "Anti-CD38 refractory" = daratumumab_or_isatuximab_status == 1,
    "Triple refractory" = triple_class_refractory_i_mi_d_pi_anti_cd38 == 1,
    "Penta refractory" = penta_refractory_2_p_is_2_i_mi_ds_anti_cd38 == 1,
    "Prior BCMA-directed therapy" = !is.na(most_recent_bcma_agent_name),
    "Most recent BCMA" = factor(case_when(
      most_recent_bcma_agent_name %in% c("Abecma", "ABECMA","Abecma on trial","Ide-cel") ~ "CAR-T",
      most_recent_bcma_agent_name %in% c("Carvykti","carvykti", "Carvikty", "Clita-cel", "Cilta-cel") ~ "CAR-T",
      most_recent_bcma_agent_name %in% c("CC98633 CAR-T","ddBCMA") ~ "CAR-T",
      most_recent_bcma_agent_name %in% c("BLENREP", "blenrep", "Belantamab") ~ "ADC",
      most_recent_bcma_agent_name %in% c("Teclistamab") ~ "TCE"
    ), levels = c("CAR-T","TCE", "ADC")),
    "Trial eligible" = eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,
    "Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
     "Baseline ferritin (ng/mL)" = as.double(baseline_ferritin_prior_to_bs_ab_initiation),
     "Baseline Hgb" = as.double(baseline_hgb),
     "Baseline LDH" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Duration of prior BCMA" = as.numeric(ymd(date_of_last_exposure_to_prior_bcma_agent) - ymd(date_of_first_exposure_to_prior_bcma_agent)),
    bs_ab_name
  ) %>%
  tbl_summary(
    by = bs_ab_name,
    missing = "no",
    #sort = list(all_categorical() ~ "frequency"),
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)",
      "Age (years)" ~ "{median} ({min} to {max})"
    )
  ) %>%
  add_p() %>%
  bold_labels() %>%
  add_variable_grouping(
    "High-risk cytogenetic abnormality at any time" = c("del(17p)","t(4;14)", "t(14;16)", "gain/amp(1q)", "amp(1q)"),
    "Refractory status" = c("Bortezomib refractory", "Lenalidomide refractory", "Pomalidomide refractory", "Anti-CD38 refractory")
  )

Characteristic	elra N = 1301	tec N = 2771	p-value2
Age (years)	71 (39 to 95)	69 (35 to 88)	0.085
Female sex	74 (57%)	128 (46%)	0.044
Race			0.14
White	102 (82%)	211 (76%)
Black	16 (13%)	54 (20%)
AAPI	3 (2.4%)	7 (2.5%)
American Indian	2 (1.6%)	0 (0%)
Other	2 (1.6%)	4 (1.4%)
ECOG			0.5
0	22 (17%)	52 (20%)
1	60 (48%)	130 (49%)
2	30 (24%)	65 (25%)
3+	14 (11%)	18 (6.8%)
LDH (U/L)	223 (177, 277)	208 (166, 287)	0.4
Hgb (g/dL)	10.05 (8.70, 11.70)	9.60 (8.30, 11.30)	0.020
Heavy chain			<0.001
IgG	75 (58%)	122 (44%)
IgA	28 (22%)	47 (17%)
IgM	2 (1.5%)	1 (0.4%)
IgD	2 (1.5%)	1 (0.4%)
None	23 (18%)	106 (38%)
True EMD	28 (22%)	36 (13%)	0.028
Oligo or non-secretory			0.3
0	116 (91%)	259 (94%)
1	12 (9.4%)	15 (5.4%)
2	0 (0%)	2 (0.7%)
≥1 HRCA	85 (69%)	164 (63%)	0.2
High-risk cytogenetic abnormality at any time
del(17p)	31 (25%)	70 (27%)	0.6
t(4;14)	17 (13%)	34 (13%)	>0.9
t(14;16)	8 (6.3%)	15 (5.8%)	0.8
gain/amp(1q)	70 (56%)	128 (49%)	0.2
amp(1q)	9 (7.9%)	18 (6.9%)	0.7
Prior LOTs	6.00 (4.00, 8.00)	6.00 (4.00, 8.00)	0.8
Prior ASCT	66 (51%)	163 (59%)	0.13
Prior CAR-T	54 (42%)	103 (37%)	0.4
Refractory status
Bortezomib refractory	87 (67%)	226 (82%)	<0.001
Lenalidomide refractory	110 (85%)	238 (86%)	0.7
Pomalidomide refractory	96 (74%)	226 (82%)	0.062
Anti-CD38 refractory	116 (89%)	267 (97%)	0.002
Triple refractory	118 (91%)	253 (92%)	0.8
Penta refractory	64 (49%)	141 (51%)	0.7
Prior BCMA-directed therapy	64 (49%)	152 (55%)	0.3
Most recent BCMA			0.012
CAR-T	54 (84%)	57 (70%)
TCE	6 (9.4%)	5 (6.2%)
ADC	4 (6.3%)	19 (23%)
Trial eligible	28 (22%)	63 (23%)	0.8
Received full dose	121 (93%)	268 (97%)	0.093
Baseline ferritin (ng/mL)	293 (132, 832)	393 (131, 1,254)	0.3
Baseline Hgb	10.05 (8.70, 11.70)	9.60 (8.30, 11.30)	0.020
Baseline LDH	223 (177, 277)	208 (166, 287)	0.4
Duration of prior BCMA	0 (0, 22)	0 (0, 126)	0.034
1 Median (Min to Max); n (%); Median (Q1, Q3)
2 Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test

4.3 Unweighted and IPTW baseline characteristics

theme_gtsummary_compact()

## Setting theme "Compact"

data1 = IPTW.data %>%
  transmute(
    "Age (years)" = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female == 1,
    "Race" = factor(case_when(
      is.na(race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5) ~ NA,
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 0 ~ "White",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 1 ~ "Black",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 2 ~ "AAPI",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 3 ~ "American Indian",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 4 ~ "Other",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 5 ~ NA
    ), levels = c("White","Black","AAPI","American Indian","Other",NA)),
    #ecog_ps_at_bs_ab_initiation_0_5,
    "ECOG" = case_when(
      ecog_ps_at_bs_ab_initiation_0_5 == 0 ~ "0",
      ecog_ps_at_bs_ab_initiation_0_5 == 1 ~ "1",
      ecog_ps_at_bs_ab_initiation_0_5 == 2 ~ "2",
      ecog_ps_at_bs_ab_initiation_0_5 >= 3 ~ "3+"
    ),
    "LDH (U/L)" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Hgb (g/dL)" = baseline_hgb,
    "Baseline ferritin (ng/mL)" = as.double(baseline_ferritin_prior_to_bs_ab_initiation),
    "Heavy chain" = factor(case_when(
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgG" ~ "IgG",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgM" ~ "IgM",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgA" ~ "IgA",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgD" ~ "IgD",
      .default = "None"
      ), levels = c("IgG", "IgA","IgM","IgD","None")),
    "True EMD" = just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,
    "≥1 HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes==1 | t_4_14_prior_to_bs_ab_0_no_1_yes == 1 | t_14_16_prior_to_bs_ab_0_no_1_yes == 1 | amp_1_q_only_0_no_1_yes == 1 | gain_amp1q21_prior_to_bs_ab_0_no_1_yes == 1,1,0),
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "Prior ASCT" = prior_auto_sct_no_0_yes_1,
    "Prior CAR-T" = prior_cart_0_no_1_yes,
    "Triple refractory" = triple_class_refractory_i_mi_d_pi_anti_cd38 == 1,
    "Penta refractory" = penta_refractory_2_p_is_2_i_mi_ds_anti_cd38 == 1,
    "Prior BCMA-directed therapy" = !is.na(most_recent_bcma_agent_name),
    "Trial eligible" = eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,
    "Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
    bs_ab_name,
    weights
  )

table_unweighted <- data1 %>% select(-weights) %>%
  tbl_summary(
    by = bs_ab_name,
    missing = "no",
    #sort = list(all_categorical() ~ "frequency"),
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)",
      "Age (years)" ~ "{median} ({min} to {max})"
    )
  ) %>%
  add_p() %>%
  bold_labels()

table_IPTW <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  ) %>%
  tbl_svysummary(
    by = bs_ab_name,
    missing = "no",
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)"),
    include = c("Age (years)", "Female sex","Race","ECOG","LDH (U/L)","Hgb (g/dL)","Baseline ferritin (ng/mL)", "Heavy chain","True EMD","≥1 HRCA","Prior LOTs","Prior ASCT","Prior CAR-T","Triple refractory","Penta refractory","Trial eligible","Prior BCMA-directed therapy", "Received full dose")
  ) %>%
  add_p() %>%
    bold_labels()

## The following warnings were returned during `add_p()`:
## ! For variable `Heavy chain` (`bs_ab_name`) and "statistic" and "p.value"
##   statistics: Chi-squared approximation may be incorrect
## ! For variable `Race` (`bs_ab_name`) and "statistic" and "p.value" statistics:
##   Chi-squared approximation may be incorrect

tbl_merge(list(table_IPTW, table_unweighted), tab_spanner = c("**IPTW**", "**Unweighted**"))

Characteristic	IPTW			Unweighted
	elra N = 1231	tec N = 2481	p-value2	elra N = 1233	tec N = 2483	p-value4
Age (years)	70 (64, 76)	69 (62, 76)	0.9	71 (39 to 95)	68 (35 to 88)	0.060
Female sex	66 (54%)	114 (46%)	0.2	68 (55%)	115 (46%)	0.11
Race			0.2			0.15
White	94 (79%)	187 (76%)		95 (81%)	187 (76%)
Black	17 (14%)	50 (20%)		16 (14%)	51 (21%)
AAPI	3 (2.5%)	6 (2.3%)		3 (2.5%)	5 (2.0%)
American Indian	2 (1.8%)	0 (0%)		2 (1.7%)	0 (0%)
Other	2 (2.0%)	4 (1.8%)		2 (1.7%)	4 (1.6%)
ECOG			>0.9			0.5
0	23 (18%)	49 (20%)		21 (17%)	50 (20%)
1	62 (50%)	119 (48%)		59 (48%)	121 (49%)
2	28 (23%)	59 (24%)		29 (24%)	60 (24%)
3+	10 (8.2%)	20 (8.2%)		14 (11%)	17 (6.9%)
LDH (U/L)	223 (178, 280)	208 (167, 290)	0.4	222 (175, 277)	208 (167, 291)	0.5
Hgb (g/dL)	9.90 (8.50, 11.40)	9.80 (8.40, 11.50)	0.9	10.00 (8.70, 11.70)	9.60 (8.25, 11.20)	0.027
Baseline ferritin (ng/mL)	346 (140, 1,250)	379 (125, 1,219)	>0.9	293 (137, 832)	393 (128, 1,254)	0.3
Heavy chain			0.001			<0.001
IgG	72 (58%)	109 (44%)		72 (59%)	107 (43%)
IgA	26 (21%)	43 (17%)		26 (21%)	42 (17%)
IgM	1 (0.7%)	1 (0.4%)		1 (0.8%)	1 (0.4%)
IgD	2 (1.5%)	0 (0%)		2 (1.6%)	0 (0%)
None	23 (19%)	95 (38%)		22 (18%)	98 (40%)
True EMD	20 (16%)	39 (16%)	>0.9	26 (21%)	34 (14%)	0.067
≥1 HRCA	81 (70%)	151 (65%)	0.3	81 (70%)	152 (65%)	0.4
Prior LOTs	6.00 (5.00, 8.00)	6.00 (4.00, 8.00)	0.4	6.00 (4.00, 8.00)	6.00 (4.00, 8.00)	>0.9
Prior ASCT	66 (54%)	141 (57%)	0.6	61 (50%)	144 (58%)	0.12
Prior CAR-T	53 (43%)	90 (36%)	0.2	50 (41%)	94 (38%)	0.6
Triple refractory	112 (91%)	224 (91%)	0.9	112 (91%)	226 (91%)	>0.9
Penta refractory	61 (49%)	123 (49%)	>0.9	61 (50%)	126 (51%)	0.8
Trial eligible	24 (19%)	55 (22%)	0.5	26 (21%)	53 (21%)	>0.9
Prior BCMA-directed therapy	67 (54%)	133 (54%)	>0.9	60 (49%)	138 (56%)	0.2
Received full dose	114 (93%)	241 (97%)	0.059	114 (93%)	241 (97%)	0.045
1 Median (Q1, Q3); n (%)
2 Design-based KruskalWallis test; Pearson’s X^2: Rao & Scott adjustment
3 Median (Min to Max); n (%); Median (Q1, Q3)
4 Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test

5 STACKED BAR GRAPH

5.1 FIGURE 1A: Toxicity, IPTW

data1 = IPTW.data %>%
  transmute(
    BsAb = as.factor(bs_ab_name),
    CRS.any = ifelse( max_crs_grade_0_5 != 0, 1, 0),
    ICANS.any = ifelse( max_icans_grade != 0,1,0),
    
    CRS.two = ifelse( max_crs_grade_0_5 == 2 | max_crs_grade_0_5 == 3 | max_crs_grade_0_5 == 4, 1, 0),
    ICANS.two = ifelse( max_icans_grade ==2 | max_icans_grade == 3 | max_icans_grade == 4,1,0),
    
    CRS.three = ifelse( max_crs_grade_0_5 == 3 | max_crs_grade_0_5 == 4, 1, 0),
    ICANS.three = ifelse( max_icans_grade == 3 | max_icans_grade == 4,1,0),
    weights
  )


weighted_table <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  )


combined_data <- rbind(
  IPTW.data %>%
    count(bs_ab_name, max_crs_grade_0_5, wt=weights) %>%
    group_by(bs_ab_name) %>%
    mutate(
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      total = sum(n[max_crs_grade_0_5 != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_crs_grade_0_5 != 0]),
      Grade = factor(max_crs_grade_0_5, levels = c(0, 4, 3, 2, 1)),
      Category = "CRS"
    ),
  IPTW.data %>%
    count(bs_ab_name, max_icans_grade, wt=weights) %>%
    group_by(bs_ab_name) %>%
    mutate(
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      total = sum(n[max_icans_grade != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_icans_grade != 0]),
      Grade = factor(max_icans_grade, levels = c(0, 4, 3, 2, 1)),
      Category = "ICANS"
    )
) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )


# Any grade toxicity
CRS.pval <- as.numeric(svychisq(~CRS.any+BsAb, weighted_table)$p.value)
ICANS.pval <- as.numeric(svychisq(~ICANS.any+BsAb, weighted_table)$p.value)

f_labels <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.pval, ICANS.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )

# G2+ toxicity
CRS.two.pval <- as.numeric(svychisq(~CRS.two+BsAb, weighted_table)$p.value)
ICANS.two.pval <- as.numeric(svychisq(~ICANS.two+BsAb, weighted_table)$p.value)

f_labels_two <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.two.pval, ICANS.two.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )

# G3+ toxicity
CRS.three.pval <- as.numeric(svychisq(~CRS.three+BsAb, weighted_table)$p.value)
ICANS.three.pval <- as.numeric(svychisq(~ICANS.three+BsAb, weighted_table)$p.value)

f_labels_three <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.three.pval, ICANS.three.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )




toxplot_IPTW <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = Grade) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(Grade != 0)) +
    geom_text(aes(label = paste0(signif(pct,2), "% (", round(n,0) ,")")), 
            data = . %>% filter(Grade != 0 & Grade != 4),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (", signif(total,2),")"), x = bs_ab_name, y = total_pct + 6, vjust = 0), 
             data = . %>% filter(Grade == 1),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_fill_brewer(palette = "Pastel1") +
  #guides(fill = guide_legend(reverse = TRUE)) +
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "right"
  ) + 
  coord_cartesian(ylim = c(0, 115)) +
  geom_text( aes(label = paste0("Grade 1+, p=", signif(label,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) + 
  geom_text( aes(label = paste0("Grade 2+, p=", signif(label,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_two ) + 
  geom_text( aes(label = paste0("Grade 3+, p=", signif(label,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_three ) + 
  ggtitle("A) IPTW")

toxplot_IPTW

5.2 FIGURE 1B: Toxicity, unweighted

combined_data <- rbind(
  combined_db %>%
    count(bs_ab_name, max_crs_grade_0_5) %>%
    group_by(bs_ab_name) %>%
    mutate(
      total = sum(n[max_crs_grade_0_5 != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_crs_grade_0_5 != 0]),
      Grade = factor(max_crs_grade_0_5, levels = c(0, 4, 3, 2, 1)),
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      Category = "CRS"
    ),
  combined_db %>%
    count(bs_ab_name, max_icans_grade) %>%
    group_by(bs_ab_name) %>%
    mutate(
      total = sum(n[max_icans_grade != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_icans_grade != 0]),
      Grade = factor(max_icans_grade, levels = c(0, 4, 3, 2, 1)),
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      Category = "ICANS"
    )
) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
    
  )




## Any-grade toxicity
CRS.d.test <- data.frame(
  results = ifelse(combined_db$max_crs_grade_0_5 !=0,1,0) ,
  test = combined_db$bs_ab_name
)
CRS.pval <- chisq.test(CRS.d.test$results,CRS.d.test$test, correct = FALSE)$p.value

ICANS.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade !=0,1,0) ,
  test = combined_db$bs_ab_name
)
ICANS.pval <- chisq.test(ICANS.d.test$results,ICANS.d.test$test, correct = FALSE)$p.value


f_labels <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.pval, ICANS.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )


## G2+ toxicity
CRS.two.d.test <- data.frame(
  results = ifelse(combined_db$max_crs_grade_0_5 == 2 | combined_db$max_crs_grade_0_5 == 3 | combined_db$max_crs_grade_0_5 == 4,1,0),
  test = combined_db$bs_ab_name
)
CRS.two.pval <- chisq.test(CRS.two.d.test$results,CRS.two.d.test$test, correct = FALSE)$p.value

ICANS.two.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade == 2 | combined_db$max_icans_grade == 3 | combined_db$max_icans_grade == 4,1,0),
  test = combined_db$bs_ab_name
)
ICANS.two.pval <- chisq.test(ICANS.two.d.test$results,ICANS.two.d.test$test, correct = FALSE)$p.value

f_labels_two <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.two.pval, ICANS.two.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )

## G3+ toxicity
CRS.three.d.test <- data.frame(
  results = ifelse( combined_db$max_crs_grade_0_5 == 3 | combined_db$max_crs_grade_0_5 == 4,1,0),
  test = combined_db$bs_ab_name
)
CRS.three.pval <- fisher.test(CRS.three.d.test$results,CRS.three.d.test$test)$p.value

ICANS.three.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade == 3 | combined_db$max_icans_grade == 4,1,0),
  test = combined_db$bs_ab_name
)
ICANS.three.pval <- chisq.test(ICANS.three.d.test$results,ICANS.three.d.test$test, correct = FALSE)$p.value

f_labels_three <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.three.pval, ICANS.three.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )



toxplot_unweighted <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = Grade) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(Grade != 0)) +
    geom_text(aes(label = paste0(signif(pct,2), "% (",n,")")), 
            data = . %>% filter(Grade != 0 & Grade != 4),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (",total,")"), x = bs_ab_name, y = total_pct + 6, vjust = 0), 
             data = . %>% filter(Grade == 1),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 100)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "none"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  scale_fill_brewer(palette = "Pastel1", direction = 1) +
  geom_text( aes(label = paste0("Grade 1+, p=", signif(label,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) + 
  geom_text( aes(label = paste0("Grade 2+, p=", signif(label,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_two ) + 
  geom_text( aes(label = paste0("Grade 3+, p=", signif(label,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_three ) + 
  ggtitle("B) Unweighted")

toxplot_unweighted

5.3 ***FIGURE 1: IPTW and unweighted toicity

FIGURE1 <- toxplot_IPTW / toxplot_unweighted

ggsave("svg/FIGURE1.svg", FIGURE1, device = "svg",
       width = 7,
       height = 11,
       units = c("in"))

ggsave("tiff/FIGURE1.tiff", FIGURE1, device = "tiff",
       width = 7,
       height = 11,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

ggsave("jpeg/FIGURE1.jpeg", FIGURE1, device = "jpeg",
       width = 7,
       height = 11,
       dpi = 600,
       units = "in")

knitr::include_graphics("jpeg/FIGURE1.jpeg")

5.4 FIGURE 2A: Response, IPTW

data1 = IPTW.data %>%
  filter(!is.na(best_response)) %>%
  transmute(
    BsAb = as.factor(bs_ab_name),
    best_response = factor(best_response, levels =  c("sCR", "CR","VGPR","PR","MR","SD","PD")),
    ORR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR",1,0),
    VGPR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR",1,0),
    CR = ifelse(best_response == "sCR" | best_response == "CR",1,0),
    weights
  )

weighted_table <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  )



data1 <-   IPTW.data %>%
  filter(!is.na(best_response)) %>%
  count(bs_ab_name, best_response, wt=weights) %>%
  group_by(bs_ab_name) %>%
  mutate(
    total = sum(n[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    pct = prop.table(n) * 100,
    total_pct = sum(pct[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    best_response = factor(best_response, levels =  c("sCR", "CR","VGPR","PR","MR", "SD","PD")),
    bs_ab_name = as.factor(bs_ab_name),
    Category = "Response rate"
  )



Response.pval <- as.numeric(svychisq(~ORR+BsAb, weighted_table)$p.value)
CR.pval <- as.numeric(svychisq(~CR+BsAb, weighted_table)$p.value)
VGPR.pval <- as.numeric(svychisq(~VGPR+BsAb, weighted_table)$p.value)

f_labels <- data.frame(Category = c("Response rate"), label = c(Response.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_CR <- data.frame(Category = c("Response rate"), label = c(CR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_VGPR <- data.frame(Category = c("Response rate"), label = c(VGPR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

responseplot_IPTW <- ggplot(data1) +
  aes(x = bs_ab_name, y = pct, fill = best_response) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR")) +
    geom_text(aes(label = paste0(signif(pct,2), "% (", round(n,0) ,")")), 
            data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR"| best_response == "PR"),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (", signif(total,2),")"), x = bs_ab_name, y = total_pct + 7, vjust = 0), 
             data = . %>% filter(best_response == "PR"),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_fill_brewer(
    palette = "Pastel1",
    direction = -1,
    breaks = rev(levels(data1$best_response))
  ) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "right"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  geom_text( aes(label = paste0("ORR, p=", signif(Response.pval,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels  ) + 
  geom_text( aes(label = paste0("≥VGPR, p=", signif(VGPR.pval,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_VGPR  ) + 
  geom_text( aes(label = paste0("≥CR, p=", signif(CR.pval,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_CR  ) + 
  labs(fill = "Best response")+ 
  ggtitle("A) IPTW")

responseplot_IPTW

5.5 FIGURE 2B: Response, unweighted

data1 <- combined_db %>%
  filter(!is.na(best_response)) %>%
  transmute(
    best_response = factor(
      best_response, 
      levels =  c("sCR", "CR","VGPR", "PR","MR","SD","PD")),
    ORR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" |best_response == "PR",1,0),
    bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
    Category = "Response rate"
  )

combined_data <- combined_db %>%
  filter(!is.na(best_response)) %>%
  count(bs_ab_name, best_response) %>%
  group_by(bs_ab_name) %>%
  transmute(
    bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
    best_response = factor(
      best_response, 
      levels =  c("sCR", "CR","VGPR", "PR","MR","SD","PD")),
    total = sum(n[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    pct = prop.table(n) * 100,
    n,
    total_pct = sum(pct[best_response == "sCR" | best_response == "CR" |best_response == "VGPR" | best_response == "PR"]),
    Category = "Response rate"
  )

## ORR
Response.d.test <- data.frame(
  results = ifelse(data1$ORR,1,0),
  test = data1$bs_ab_name
)
Response.pval <- chisq.test(Response.d.test$results,Response.d.test$test, correct = FALSE)$p.value

## VGPR
VGPR.d.test <- data.frame(
  results = ifelse(data1$best_response == "CR" | data1$best_response == "VGPR",1,0),
  test = data1$bs_ab_name
)
VGPR.pval <- chisq.test(VGPR.d.test$results,VGPR.d.test$test, correct = FALSE)$p.value

## CR
CR.d.test <- data.frame(
  results = ifelse(data1$best_response == "CR",1,0),
  test = data1$bs_ab_name
)
CR.pval <- chisq.test(CR.d.test$results,CR.d.test$test, correct = FALSE)$p.value


f_labels <- data.frame(Category = c("Response rate"), label = c(Response.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_CR <- data.frame(Category = c("Response rate"), label = c(CR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_VGPR <- data.frame(Category = c("Response rate"), label = c(VGPR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )


responseplot_unweighted <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = best_response) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(best_response == "sCR" | best_response == "CR" |best_response == "VGPR" | best_response == "PR")) +
  geom_text(aes(label = paste0(signif(pct,2), "% (",n,")")), 
             data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"),
             position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (",total,")"), x = bs_ab_name, y = total_pct + 7, vjust = 0), 
             data = . %>% filter(best_response == "PR"),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "none"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  scale_fill_brewer(palette = "Pastel1", direction = -1) +
  geom_text( aes(label = paste0("ORR, p=", signif(Response.pval,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels) + 
  geom_text( aes(label = paste0("≥VGPR, p=", signif(VGPR.pval,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_VGPR  ) + 
  geom_text( aes(label = paste0("≥CR, p=", signif(CR.pval,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_CR  ) + 
  labs(fill = "Best response")+ 
  ggtitle("B) Unweighted")

responseplot_unweighted

5.6 ***FIGURE 2: IPTW and unweighted response

FIGURE2 <- responseplot_IPTW / responseplot_unweighted

ggsave("svg/FIGURE2.svg", FIGURE2, device = "svg",
       width = 7,
       height = 11,
       units = c("in"))

ggsave("tiff/FIGURE2.tiff", FIGURE2, device = "tiff",
       width = 7,
       height = 11,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

ggsave("jpeg/FIGURE2.jpeg", FIGURE2, device = "jpeg",
       width = 7,
       height = 11,
       dpi = 600,
       units = "in")

knitr::include_graphics("jpeg/FIGURE2.jpeg")

6 SURVIVAL ANALYSIS

6.1 Length of follow-up (Reverse KM, months)

data1 <- IPTW.data %>%
  transmute(
    bs_ab_name,
    Death = ifelse(death_yes_no,1,0),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Reverse_death = ifelse(death_yes_no == 1, 0,1),
    weights
  )

reverse_km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1)
reverse_km_OS

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1, weights = weights)
## 
##    15 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     356 356    215   13.4    12.5    14.9

reverse_km_OS_elra <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1 %>% filter(bs_ab_name == "elra"))
reverse_km_OS_elra

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1 %>% filter(bs_ab_name == "elra"), weights = weights)
## 
##    6 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     117 117   74.5   7.47    6.03       9

reverse_km_OS_tec <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1 %>% filter(bs_ab_name == "tec"))
reverse_km_OS_tec

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1 %>% filter(bs_ab_name == "tec"), weights = weights)
## 
##    9 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     239 239    140   17.1    14.9    20.1

6.2 Stratified by BsAb

data1 <- IPTW.data %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

data_dor <- IPTW.data %>%
  filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
  transmute(
    weights,
    bs_ab_name,
    day_30_response,
    day_90_response,
    best_response,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_of_bs_ab_step_up_d1,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    day_30_response,
    day_90_response,
    best_response,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
    ),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
  ) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)

km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, 
##     data = data1, weights = weights)
## 
##    15 observations deleted due to missingness 
##                 records   n events median 0.95LCL 0.95UCL
## bs_ab_name=elra     117 117   42.8     NA     8.6      NA
## bs_ab_name=tec      239 239   98.6   21.5    16.4      NA

med_time <- surv_median(km_OS)

cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))

OS <- ggsurvplot(km_OS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Overall survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "A)",
           size = 0.5
) 

OS$plot <- OS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

OS

km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS

## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data1, weights = weights)
## 
##    12 observations deleted due to missingness 
##                 records   n events median 0.95LCL 0.95UCL
## bs_ab_name=elra     118 118   68.1   3.77    2.57    7.17
## bs_ab_name=tec      241 241  151.1   8.50    7.10   11.30

med_time <- surv_median(km_PFS)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))


PFS <- ggsurvplot(km_PFS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Progression-free survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "B)",
           size = 0.5
)

PFS$plot <- PFS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

PFS

km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR

## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data_dor, weights = weights)
## 
##    3 observations deleted due to missingness 
##                 records     n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      62  60.2   20.6   11.7    6.13      NA
## bs_ab_name=tec      118 120.7   56.1   13.7    9.27    22.3

med_time <- surv_median(km_DOR)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))


DOR <- ggsurvplot(km_DOR,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           #risk.table.fontsize = 4,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Duration of response",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "C)",
           size = 0.5
)

DOR$plot <- DOR$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

DOR

6.3 ***FIGURE 3: PFS/OS/DOR, stratified by trial elibility

FIGURE3 <- wrap_plots(
  (OS$plot / OS$table) + plot_layout(heights = c(4, 0.9)),
  (PFS$plot / PFS$table) + plot_layout(heights = c(4, 0.9)),
  (DOR$plot / DOR$table) + plot_layout(heights = c(4, 0.9)),
  ncol = 1  # Force vertical stacking
)

ggsave("svg/FIGURE3.svg", FIGURE3, device = "svg",
       width = 6.5,
       height = 8.5,
       units = c("in"))

ggsave("jpeg/FIGURE3.jpeg", FIGURE3, device = "jpeg",
       width = 6.5,
       height = 8.5,
       dpi = 600,
       units = c("in"))

ggsave("tiff/FIGURE3.tiff", FIGURE3, device = "tiff",
       width = 6.5,
       height = 8.5,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

knitr::include_graphics("jpeg/FIGURE3.jpeg")

6.4 Schoenfeld residuals and PH test

data1 <- IPTW.data %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    elra_ind = ifelse(bs_ab_name == "elra", 1, 0),  # 1 = elra, 0 = tec
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

## ---- Cox models ----
cox_OS <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ elra_ind, weights = weights, data = data1)

cox_PFS <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ elra_ind, weights = weights, data = data1)

## ---- PH diagnostics ----
zph_OS  <- cox.zph(cox_OS)
zph_PFS <- cox.zph(cox_PFS)

# Global PH p-values (same as elra_ind right now)
p_OS  <- signif(zph_OS$table["GLOBAL",  "p"], 3)
p_PFS <- signif(zph_PFS$table["GLOBAL", "p"], 3)

par(mfrow = c(1, 2))

plot(
  zph_OS,
  main = paste0(
    "Schoenfeld residuals – OS\n",
    "Global PH test: p = ", p_OS
  )
)

plot(
  zph_PFS,
  main = paste0(
    "Schoenfeld residuals – PFS\n",
    "Global PH test: p = ", p_PFS
  )
)

par(mfrow = c(1, 1))

7 REGRESSION

7.1 ORR, ≥CR, ≥VGPR, PFS, OS, DOR, IFS, CRS G1+, ICANS G1+, steroid use (IPTW)

options(gtsummary.add_p_footnote = FALSE, gtsummary.add_estimate_to_reference_rows = FALSE)
theme_gtsummary_compact()

## Setting theme "Compact"

preds0 <- c("Response", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds1 <- c("CR/sCR", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds2 <- c("BsAb", "Baseline Hgb",  "Baseline LDH", "ECOG 2+","True EMD", "Prior LOTs", "Prior BCMA")
preds3 <- c("ICANS 1+", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds4 <- c("CRS 1+", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds5 <- c("Steroid use", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")

preds6 <- c("VGPR", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")

data0 <- IPTW.data %>%
  transmute(
    weights,
    BsAb = factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra"), levels = c("Tec","Elra")),
    HT_Risk,
    best_response,
    "Infection" = ifelse(infection == 1, 1,0),
    "CrCl <40" = ifelse(cr_cl_prior_to_step_up_dose_number_1 <40, 1,0),
    "CR/sCR" = ifelse(best_response %in% c("CR","sCR"),1,0),
    "VGPR" = ifelse(best_response %in% c("CR","sCR", "VGPR"),1,0),
    "Response" = ifelse(best_response %in% c("CR","sCR", "VGPR", "PR"),1,0),
    "True EMD" = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2, 1,0),
    "Steroid use" = ifelse(steroid_use_as_tx==1,1,0),
    "Age" = age_at_bs_ab_initiation/10,
    "ECOG" = ecog_ps_at_bs_ab_initiation_0_5,
    "ECOG 2+" = case_when(
      is.na(ecog_ps_at_bs_ab_initiation_0_5) ~ NA,
      ecog_ps_at_bs_ab_initiation_0_5 %in% c(2,3,4) ~ 1,
      TRUE ~ 0),
    "ICANS 1+" = case_when(
      is.na(max_icans_grade) ~ NA,
      max_icans_grade %in% c(1,2,3,4) ~1,
      TRUE ~ 0
      ),
    "CRS 1+" = case_when(
      is.na(max_crs_grade_0_5) ~ NA,
      max_crs_grade_0_5 %in% c(1,2,3,4) ~1,
      TRUE ~ 0
      ),
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes | t_4_14_prior_to_bs_ab_0_no_1_yes | t_14_16_prior_to_bs_ab_0_no_1_yes | gain_amp1q21_prior_to_bs_ab_0_no_1_yes, 1, 0),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no,
    prior_bcma_directed_therapy_yes_no,
    Days.since.last.BCMA = case_when(
      is.na(prior_bcma_directed_therapy_yes_no) ~ NA,
      date_of_last_exposure_to_prior_bcma_agent > date_of_bs_ab_step_up_d1 ~ NA,
      prior_bcma_directed_therapy_yes_no == 1 ~ as.numeric(as.Date(date_of_bs_ab_step_up_d1) - as.Date(date_of_last_exposure_to_prior_bcma_agent))
    ),
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    "ALC BL" = baseline_alc,
    "CrCl BL" = cr_cl_prior_to_step_up_dose_number_1,
    "Trial ineligible" = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,0,1),
    "Baseline platelets" = baseline_pl_ts,
    "Baseline Hgb" = baseline_hgb,
    "Baseline LDH" = as.double(baseline_ldh_prior_to_bs_ab_initiation)/100,
    "Baseline ferritin" = as.double(baseline_ldh_prior_to_bs_ab_initiation)/1000,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
      ),
    date_of_bs_ab_step_up_d1,
    bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes,
    days_to_infection_or_last_contact = as.numeric(
      pmin(as.Date(date_of_infection_74), as.Date(date_of_last_contact), na.rm = TRUE) - as.Date(date_of_bs_ab_step_up_d1)
    ),
    infection_or_death = ifelse(!is.na(date_of_infection_74) | death_yes_no == 1,1,0),
    IVIg = ifelse(ivig_use_yes_1_no_0 == 1, 1, 0),
    ALPS = ifelse(baseline_hgb >10 & baseline_ldh_prior_to_bs_ab_initiation <300,1,0),
    baseline_pl_ts,
    baseline_anc,
    baseline_hgb,
    baseline_crp_prior_to_bs_ab_initiation,
    baseline_ferritin_prior_to_bs_ab_initiation,
    CAR.HEMATOTOX = baseline_pl_ts * -0.05 + baseline_anc * -0.03 + baseline_hgb * -0.20 + as.numeric(baseline_crp_prior_to_bs_ab_initiation) * 0.15 + as.numeric(baseline_ferritin_prior_to_bs_ab_initiation) * 0.001
    ) %>%
  mutate(
    "BCMA exposure status" = factor(case_when(
      #prior_bcma_directed_therapy_yes_no == 0 ~ "  Never",
      Days.since.last.BCMA < 365 ~ "  <1 yr",
      Days.since.last.BCMA >= 365 ~ "  ≥1 yr"
    ), c("  ≥1 yr","  <1 yr")),
    `CAR-HEMATOTOX risk` = factor(
      case_when(
        CAR.HEMATOTOX < 0 ~ "Low",
        CAR.HEMATOTOX <= 5 ~ "Intermediate",
        CAR.HEMATOTOX > 5 ~ "High"
      ),
      levels = c("Low", "Intermediate", "High")
  )
  )

data1 <- data0 %>% filter(bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes==0) ## Exclude patients who received bsAb as bridging

data2 <- data1 %>% 
  filter(best_response %in% c("PR","VGPR","sCR", "CR")) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)


## Response
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds0)),
  method = glm,
  y = `Response`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")

## There was a warning constructing the model for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "ECOG 2+". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "True EMD". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "ECOG 2+". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "True EMD". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!

mv_tab<-glm( `Response` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds0), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**") 

table_response <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## VGPR or better
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds6)),
  method = glm,
  y = `VGPR`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE,
) %>% modify_header(estimate = "**Estimate**")

## There was a warning constructing the model for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "ECOG 2+". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "True EMD". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "ECOG 2+". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "True EMD". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!

mv_tab<-glm( `VGPR` ~ BsAb +`True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds6), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_VGPR <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## sCR/CR
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds1)),
  method = glm,
  y = `CR/sCR`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE,
) %>% modify_header(estimate = "**Estimate**")

## There was a warning constructing the model for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "ECOG 2+". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "True EMD". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "ECOG 2+". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "True EMD". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!

mv_tab<-glm( `CR/sCR` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds1), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_sCR.CR <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## OS

uv_tab <- tbl_uvregression(
  data1[c(preds2)],
  method = coxph,
  y = Surv(data1$Days.to.death.or.DLC, data1$Death),
  exponentiate = TRUE,
  method.args = list(weights = data1$weights)
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(Days.to.death.or.DLC, Death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data1,  weights = weights) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_os <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 


## PFS
uv_tab <- tbl_uvregression(
  data1[c(preds2)],
  method = coxph,
  y = Surv(data1$Days.to.death.or.relapse.or.DLC, data1$Relapse.or.death),
  exponentiate = TRUE,
  method.args = list(weights = data1$weights)
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(Days.to.death.or.relapse.or.DLC, Relapse.or.death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data1) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_pfs <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
)

## DOR
uv_tab <- tbl_uvregression(
  data2[c(preds2)],
  method = coxph,
  y = Surv(data2$Days.to.death.or.relapse.or.DLC, data2$Relapse.or.death),
  exponentiate = TRUE,
  method.args = list(weights = data2$weights)
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(  Days.to.death.or.relapse.or.DLC, Relapse.or.death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data2) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_dor <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
)

## IFS
uv_tab <- tbl_uvregression(
  data1[c(preds2)],
  method = coxph,
  y = Surv(data1$days_to_infection_or_last_contact, data1$infection_or_death),
  exponentiate = TRUE,
  method.args = list(weights = data1$weights)
 ) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(  days_to_infection_or_last_contact, infection_or_death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data1) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_ifs <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
)

## ICANS 1+
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds3)),
  method = glm,
  y = `ICANS 1+`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")

## There was a warning constructing the model for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "ECOG 2+". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "True EMD". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "ECOG 2+". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "True EMD". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!

mv_tab<-glm( `ICANS 1+` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds3), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_ICANS <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## CRS 1+
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds4)),
  method = glm,
  y = `CRS 1+`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")

## There was a warning constructing the model for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "ECOG 2+". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "True EMD". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "ECOG 2+". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "True EMD". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!

mv_tab<-glm( `CRS 1+` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds4), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_CRS <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## Steroid use
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds5)),
  method = glm,
  y = `Steroid use`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")

## There was a warning constructing the model for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "ECOG 2+". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "True EMD". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "ECOG 2+". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "True EMD". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!

mv_tab<-glm( `Steroid use` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds5), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_steroid <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 


# Combine all tables into one
final_table <- tbl_stack(
  list(
    table_response %>% modify_header() %>% modify_footnote(everything() ~ NA, abbreviation = TRUE) %>%
      modify_footnote(estimate_1 ~ "Estimate = Odds ratio or hazard ratio", abbreviation = TRUE),
    table_VGPR,
    table_sCR.CR,
    table_os,
    table_pfs,
    table_dor,
    table_ifs,
    table_ICANS,
    table_CRS,
    table_steroid
  ),
  group_header = c(
    "Response",
    "≥VGPR",
    "CR/sCR",
    "Overall Survival",
    "Progression-free Survival",
    "Duration of Response",
    "Infection-free Survival",
    "ICANS 1+",
    "CRS 1+",
    "Steroid use"
  )
)
  
final_table

Characteristic	Univariate				Multivariable
	N	Estimate1	95% CI	p-value	Estimate	95% CI	p-value
Response
BsAb	370
Tec		—	—		—	—
Elra		0.86	0.55, 1.35	0.5	0.86	0.52, 1.42	0.6
Baseline Hgb	370	1.46	1.29, 1.65	<0.001	1.34	1.18, 1.53	<0.001
Baseline LDH	370	0.86	0.77, 0.94	0.002	0.86	0.77, 0.95	0.004
ECOG 2+	370	0.39	0.25, 0.61	<0.001	0.42	0.26, 0.69	<0.001
True EMD	370	0.91	0.52, 1.63	0.8	0.94	0.50, 1.78	0.9
Prior LOTs	370	0.87	0.79, 0.94	0.001	0.91	0.82, 1.01	0.087
Prior BCMA	370	0.56	0.36, 0.85	0.008	0.50	0.29, 0.87	0.014
≥VGPR
BsAb	370
Tec		—	—		—	—
Elra		0.62	0.40, 0.95	0.031	0.57	0.35, 0.93	0.026
Baseline Hgb	370	1.38	1.23, 1.55	<0.001	1.32	1.17, 1.49	<0.001
Baseline LDH	370	0.84	0.73, 0.93	0.003	0.84	0.73, 0.93	0.004
ECOG 2+	370	0.52	0.33, 0.81	0.004	0.57	0.34, 0.93	0.026
True EMD	370	0.63	0.35, 1.10	0.11	0.56	0.30, 1.04	0.071
Prior LOTs	370	0.88	0.80, 0.95	0.003	0.95	0.86, 1.06	0.4
Prior BCMA	370	0.47	0.31, 0.72	<0.001	0.41	0.24, 0.70	0.001
CR/sCR
BsAb	370
Tec		—	—		—	—
Elra		0.70	0.44, 1.10	0.12	0.68	0.40, 1.12	0.13
Baseline Hgb	370	1.44	1.29, 1.63	<0.001	1.39	1.23, 1.58	<0.001
Baseline LDH	370	0.86	0.74, 0.96	0.019	0.87	0.75, 0.98	0.044
ECOG 2+	370	0.51	0.31, 0.81	0.006	0.61	0.36, 1.03	0.068
True EMD	370	0.66	0.35, 1.20	0.2	0.62	0.31, 1.19	0.2
Prior LOTs	370	0.89	0.81, 0.97	0.014	0.95	0.85, 1.06	0.4
Prior BCMA	370	0.56	0.37, 0.86	0.009	0.52	0.31, 0.89	0.018
Overall Survival
BsAb	352
Tec		—	—		—	—
Elra		1.41	0.96, 2.07	0.083	1.46	0.97, 2.20	0.073
Baseline Hgb	352	0.73	0.66, 0.81	<0.001	0.77	0.69, 0.86	<0.001
Baseline LDH	352	1.13	1.06, 1.20	<0.001	1.15	1.10, 1.20	<0.001
ECOG 2+	352	2.25	1.61, 3.16	<0.001	1.92	1.32, 2.78	<0.001
True EMD	352	1.34	0.88, 2.03	0.2	1.28	0.82, 2.00	0.3
Prior LOTs	352	1.06	0.99, 1.13	0.091	1.04	0.97, 1.11	0.3
Prior BCMA	352	1.09	0.77, 1.54	0.6	1.20	0.81, 1.78	0.4
Progression-free Survival
BsAb	355
Tec		—	—		—	—
Elra		1.43	1.05, 1.94	0.021	1.52	1.11, 2.09	0.009
Baseline Hgb	355	0.77	0.71, 0.83	<0.001	0.81	0.74, 0.89	<0.001
Baseline LDH	355	1.18	1.08, 1.28	<0.001	1.17	1.08, 1.27	<0.001
ECOG 2+	355	1.63	1.23, 2.16	<0.001	1.43	1.04, 1.95	0.027
True EMD	355	1.04	0.74, 1.47	0.8	0.94	0.67, 1.32	0.7
Prior LOTs	355	1.10	1.05, 1.16	<0.001	1.06	1.00, 1.12	0.056
Prior BCMA	355	1.58	1.18, 2.10	0.002	1.63	1.17, 2.27	0.004
Duration of Response
BsAb	177
Tec		—	—		—	—
Elra		1.49	0.87, 2.56	0.14	1.50	0.85, 2.63	0.2
Baseline Hgb	177	0.84	0.74, 0.95	0.007	0.83	0.72, 0.96	0.014
Baseline LDH	177	1.05	0.90, 1.23	0.5	1.04	0.89, 1.21	0.7
ECOG 2+	177	1.05	0.64, 1.72	0.8	0.89	0.52, 1.52	0.7
True EMD	177	1.48	0.75, 2.93	0.3	1.36	0.68, 2.72	0.4
Prior LOTs	177	1.01	0.89, 1.15	0.8	1.02	0.90, 1.17	0.7
Prior BCMA	177	1.08	0.67, 1.74	0.8	1.06	0.65, 1.75	0.8
Infection-free Survival
BsAb	358
Tec		—	—		—	—
Elra		1.25	0.93, 1.67	0.13	1.24	0.92, 1.67	0.2
Baseline Hgb	358	0.87	0.81, 0.94	<0.001	0.90	0.83, 0.97	0.010
Baseline LDH	358	1.10	1.05, 1.15	<0.001	1.09	1.05, 1.14	<0.001
ECOG 2+	358	1.70	1.30, 2.22	<0.001	1.52	1.14, 2.02	0.004
True EMD	358	1.16	0.83, 1.62	0.4	1.10	0.79, 1.54	0.6
Prior LOTs	358	0.98	0.93, 1.03	0.4	0.98	0.93, 1.04	0.5
Prior BCMA	358	0.85	0.65, 1.11	0.2	0.99	0.73, 1.35	>0.9
ICANS 1+
BsAb	370
Tec		—	—		—	—
Elra		0.87	0.49, 1.52	0.6	0.91	0.49, 1.65	0.8
Baseline Hgb	370	0.83	0.72, 0.95	0.008	0.84	0.72, 0.98	0.024
Baseline LDH	370	1.05	0.95, 1.15	0.3	1.04	0.93, 1.14	0.5
ECOG 2+	370	3.03	1.77, 5.21	<0.001	2.44	1.39, 4.31	0.002
True EMD	370	1.40	0.69, 2.69	0.3	1.31	0.62, 2.62	0.5
Prior LOTs	370	0.88	0.77, 0.98	0.034	0.88	0.76, 1.01	0.069
Prior BCMA	370	0.61	0.36, 1.03	0.067	1.01	0.53, 1.92	>0.9
CRS 1+
BsAb	370
Tec		—	—		—	—
Elra		0.53	0.34, 0.82	0.004	0.55	0.35, 0.86	0.010
Baseline Hgb	370	0.94	0.85, 1.04	0.3	0.91	0.81, 1.01	0.076
Baseline LDH	370	0.98	0.90, 1.06	0.5	0.97	0.89, 1.06	0.5
ECOG 2+	370	1.08	0.70, 1.68	0.7	0.97	0.60, 1.54	0.9
True EMD	370	0.58	0.33, 1.02	0.061	0.53	0.29, 0.94	0.033
Prior LOTs	370	0.92	0.84, 1.00	0.045	0.87	0.79, 0.96	0.006
Prior BCMA	370	1.17	0.77, 1.76	0.5	1.84	1.12, 3.04	0.016
Steroid use
BsAb	355
Tec		—	—		—	—
Elra		0.83	0.48, 1.40	0.5	0.87	0.50, 1.50	0.6
Baseline Hgb	355	0.89	0.79, 1.01	0.074	0.90	0.79, 1.03	0.14
Baseline LDH	355	1.07	0.98, 1.17	0.12	1.05	0.96, 1.15	0.3
ECOG 2+	355	1.71	1.03, 2.83	0.036	1.32	0.78, 2.23	0.3
True EMD	355	1.61	0.85, 2.97	0.13	1.57	0.81, 2.94	0.2
Prior LOTs	355	0.88	0.79, 0.98	0.028	0.94	0.83, 1.06	0.3
Prior BCMA	355	0.47	0.29, 0.78	0.003	0.65	0.36, 1.16	0.15
1 Estimate = Odds ratio or hazard ratio

8 CONTOUR PLOTS

8.1 FIGURE 4: Baseline Hgb and LDH

data1 <- IPTW.data %>%
  filter(bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes==0) %>% ## Exclude patients who received bsAb as bridging
  #filter(bs_ab_name == "tec") %>%
  transmute(
    bs_ab_name,
    weights,
    "CR/sCR" = ifelse(best_response %in% c("CR", "sCR"),1,0),
    "PD" = ifelse(best_response %in% c("PD"),1,0),
    "True EMD" = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2, 1,0),
    "Age" = age_at_bs_ab_initiation/10,
    "ECOG" = ecog_ps_at_bs_ab_initiation_0_5,
    ECOG_2_4 = case_when(
      is.na(ecog_ps_at_bs_ab_initiation_0_5) ~ NA,
      ecog_ps_at_bs_ab_initiation_0_5 %in% c(2,3,4) ~ 1,
      TRUE ~ 0),
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes | t_4_14_prior_to_bs_ab_0_no_1_yes | t_14_16_prior_to_bs_ab_0_no_1_yes | gain_amp1q21_prior_to_bs_ab_0_no_1_yes, 1, 0),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    "ALC BL" = baseline_alc,
    "CrCl BL" = cr_cl_prior_to_step_up_dose_number_1,
    "Trial ineligible" = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,0,1),
    "Baseline platelets" = baseline_pl_ts,
    baseline_hgb,
    baseline_ldh = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    baseline_ferritin = as.numeric(baseline_ferritin_prior_to_bs_ab_initiation),
    ) %>% filter(!is.na(baseline_ldh))


## Tec
data_tec = data1 %>% filter(bs_ab_name == "tec")

dd_tec <- datadist(data_tec)
options(datadist='dd_tec')

S_tec <- Surv(data_tec$Days.to.death.or.relapse.or.DLC/30, data_tec$Relapse.or.death)
f_tec <- cph(S_tec ~ rcs(baseline_hgb, 3) * rcs(baseline_ldh, 3) , x=TRUE, y=TRUE,surv=TRUE,data=data_tec, weights=weights)
med_tec <- Quantile(f_tec)

pred_tec  <- data.frame(Predict(f_tec,  baseline_hgb, baseline_ldh,
                                fun = function(x) med_tec(lp = x))) %>%
  mutate(drug = "Teclistamab")

## Elra
data_elra = data1 %>% filter(bs_ab_name == "elra")

dd_elra <- datadist(data_elra)
options(datadist='dd_elra')

S_elra <- Surv(data_elra$Days.to.death.or.relapse.or.DLC/30, data_elra$Relapse.or.death)
f_elra <- cph(S_elra ~ rcs(baseline_hgb, 3) * rcs(baseline_ldh, 3) , x=TRUE, y=TRUE,surv=TRUE,data=data_elra, weights=weights)
med_elra <- Quantile(f_elra)

pred_elra <- data.frame(Predict(f_elra, baseline_hgb, baseline_ldh,
                                fun = function(x) med_elra(lp = x))) %>%
  mutate(drug = "Elranatamab")

## Combine both predictions
all_pred <- bind_rows(pred_tec, pred_elra)
z_range <- range(all_pred$yhat, na.rm = TRUE)
brks    <- pretty(z_range, n = 10)   # tweak n if you want more/fewer bands

all_pred <- all_pred %>%
  mutate(
    z_band = cut(
      yhat,
      breaks        = brks,
      include.lowest = TRUE,
      right          = FALSE
    )
  )

z_levels <- levels(all_pred$z_band)

cols <- viridis_pal()(length(z_levels))


contour_plot_tec <- all_pred %>%
  filter(drug == "Teclistamab") %>%
  ggplot(aes(baseline_hgb, baseline_ldh, fill = z_band)) +
  geom_raster() +   # or geom_tile(); blocky but guarantees alignment
  theme_classic() +
  theme(
    panel.grid.major = element_line(color = "grey85"),
    panel.grid.minor = element_line(color = "grey92")
  )+
  labs(
    title = "A) Teclistamab",
    x = "Baseline Hgb (g/dL)",
    y = "Baseline LDH (U/L)",
    fill = "Median PFS (months)"
  ) +
  scale_fill_manual(
    values       = cols,
    breaks       = z_levels,
    drop         = FALSE,
    na.translate = FALSE
  ) +
  scale_x_continuous(breaks = seq(0, 20, 1))+
  scale_y_continuous(breaks = seq(0, 2000, 100))+
  coord_cartesian(xlim = c(7, 13), ylim = c(120, 780))

contour_plot_tec

contour_plot_elra <- all_pred %>%
  filter(drug == "Elranatamab") %>%
  ggplot(aes(baseline_hgb, baseline_ldh, fill = z_band)) +
  geom_raster() +
  theme_classic() +
  theme(
    legend.position = "none",
    panel.grid.major = element_line(color = "grey85"),
    panel.grid.minor = element_line(color = "grey92")
    ) +
  labs(
    title = "B) Elranatamab",
    x = "Baseline Hgb (g/dL)",
    y = "Baseline LDH (U/L)"
    #fill = "Median PFS (months)"
  ) +
  scale_fill_manual(
    values       = cols,
    breaks       = z_levels,
    drop         = FALSE,
    na.translate = FALSE
  ) +
  scale_x_continuous(breaks = seq(0, 20, 1))+
  scale_y_continuous(breaks = seq(0, 2000, 100))+
  coord_cartesian(xlim = c(7, 13), ylim = c(120, 780))

 contour_plot_elra

8.2 ***FIGURE 4: Hgb and LDH contour

combined_contour <- (contour_plot_tec / contour_plot_elra) +
  plot_layout(guides = "collect", axes = "collect")


ggsave("svg/FIGURE4.svg", combined_contour, device = "svg",
       width = 6.5,
       height = 6.5,
       units = c("in"))

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_raster()`).

## Warning: Removed 4293 rows containing missing values or values outside the scale range
## (`geom_raster()`).

ggsave("jpeg/FIGURE4.jpeg", combined_contour, device = "jpeg",
       width = 6.5,
       height = 6.5,
       dpi = 600,
       units = c("in"))

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_raster()`).
## Removed 4293 rows containing missing values or values outside the scale range
## (`geom_raster()`).

ggsave("tiff/FIGURE4.tiff", combined_contour, device = "tiff",
       width = 6.5,
       height = 6.5,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_raster()`).
## Removed 4293 rows containing missing values or values outside the scale range
## (`geom_raster()`).

knitr::include_graphics("jpeg/FIGURE4.jpeg")